In [ ]:
import datetime
# Third-party
from astropy.io import ascii
import astropy.coordinates as coord
import astropy.units as u
import astropy.time as at
import matplotlib as mpl
import matplotlib.pyplot as pl
import numpy as np
%matplotlib inline
In [ ]:
import astroplan
from astroplan import Observer, FixedTarget
from astropy.time import Time
In [ ]:
mdm = Observer.at_site("MDM", timezone="America/Phoenix")
t1 = Time(datetime.datetime(2016, 2, 15, 18, 0, tzinfo=mdm.timezone))
t2 = t1 + 12*u.hour
time_range = Time([t1, t2])
In [ ]:
def coords_in_rect(c, corner_c):
if not c.frame.is_equivalent_frame(corner_c[0].frame):
raise ValueError("Frame mismatch.")
min_lon = corner_c[0].spherical.lon
min_lat = corner_c[0].spherical.lat
max_lon = corner_c[1].spherical.lon
max_lat = corner_c[1].spherical.lat
return ((c.spherical.lon > min_lon) & (c.spherical.lon < max_lon) &
(c.spherical.lat > min_lat) & (c.spherical.lat < max_lat))
In [ ]:
css = ascii.read("/Users/adrian/projects/triand-rrlyrae/data/catalina.csv")
print(css.colnames)
In [ ]:
linear = ascii.read("/Users/adrian/projects/triand-rrlyrae/data/linear.csv")
print(linear.colnames)
In [ ]:
css_c = coord.SkyCoord(ra=css['ra']*u.deg, dec=css['dec']*u.deg, distance=css['helio_dist']*u.kpc)
lin_c = coord.SkyCoord(ra=linear['ra']*u.deg, dec=linear['dec']*u.deg, distance=linear['helio_dist']*u.kpc)
In [ ]:
fig = pl.figure(figsize=(8,6))
ax = fig.add_subplot(111)
ax.plot(css_c.galactic.l.degree, css_c.galactic.b.degree, marker='.', ls='none', alpha=0.4)
ax.plot(lin_c.galactic.l.degree, lin_c.galactic.b.degree, marker='.', ls='none', alpha=0.4)
# GASS
r = pl.Rectangle((100,15), 160, 15, zorder=-100, alpha=0.2, color='r')
ax.add_patch(r)
# A13
r = pl.Rectangle((130,20), 50, 20, zorder=-100, alpha=0.2, color='g')
ax.add_patch(r)
ax.set_xlim(360, 0)
Window and distance range taken from Rocha-pinto (2003): http://iopscience.iop.org/article/10.1086/378668/pdf
In [ ]:
window_corners = [coord.SkyCoord(l=100*u.deg, b=15*u.deg,frame='galactic'),
coord.SkyCoord(l=260*u.deg, b=35*u.deg,frame='galactic')]
css_sky_window_ix = coords_in_rect(css_c.galactic, window_corners)
lin_sky_window_ix = coords_in_rect(lin_c.galactic, window_corners)
print("{} CSS targets, {} LINEAR targets in this window.".format(css_sky_window_ix.sum(), lin_sky_window_ix.sum()))
In [ ]:
css_distance_ix = ((css_c.distance > 7*u.kpc) & (css_c.distance < 15.*u.kpc))
In [ ]:
gass_tbl = css[css_sky_window_ix & css_distance_ix]
gass = coord.SkyCoord(l=gass_tbl['l']*u.deg, b=gass_tbl['b']*u.deg,
distance=gass_tbl['helio_dist']*u.kpc, frame='galactic')
In [ ]:
print("{} GASS targets".format(len(gass)))
In [ ]:
pl.figure(figsize=(10,6))
pl.scatter(gass.l.degree, gass.b.degree, c=gass_tbl['VmagAvg'], marker='o')
pl.plot(css_c.galactic.l.degree, css_c.galactic.b.degree, marker='.', ls='none')
pl.plot(lin_c.galactic.l.degree, lin_c.galactic.b.degree, marker='.', ls='none')
pl.xlim(260,100)
pl.ylim(0,30)
In [ ]:
fig,axes = pl.subplots(1,2,figsize=(12,5))
n,bins,pa = axes[0].hist(gass.distance, bins=np.linspace(5,20,15))
axes[0].set_xlabel("Radial dist. [kpc]")
axes[0].set_xlim(bins.min(), bins.max())
n,bins,pa = axes[1].hist(gass.galactocentric.represent_as(coord.CylindricalRepresentation).rho, bins=bins+8)
axes[1].set_xlabel("Cylindrical dist. [kpc]")
axes[1].set_xlim(bins.min(), bins.max())
axes[0].set_ylabel("Number RRL")
In [ ]:
fig,ax = pl.subplots(1,1,figsize=(6,6),subplot_kw =dict(polar=True))
ax.add_artist(mpl.patches.Circle((-8.,0), radius=2.5, transform=ax.transData._b, facecolor='r', alpha=0.2))
ax.add_artist(mpl.patches.Circle((-8.,0), radius=0.5, transform=ax.transData._b, facecolor='y', alpha=1.))
gass_cyl = gass.galactocentric.represent_as(coord.CylindricalRepresentation)
css_cyl = css_c.galactocentric.represent_as(coord.CylindricalRepresentation)
ax.plot(css_cyl.phi.to(u.radian)[np.abs(css_cyl.z) < 5*u.kpc],
css_cyl.rho.to(u.kpc)[np.abs(css_cyl.z) < 5*u.kpc],
color='k', linestyle='none', marker='.', alpha=0.1)
# ax.plot(gass_cyl.phi.to(u.radian), gass_cyl.rho.to(u.kpc),
# color='k', linestyle='none', marker='o')
ax.scatter(gass_cyl.phi.to(u.radian), gass_cyl.rho.to(u.kpc),
c=gass_tbl['VmagAvg'], marker='o')
ax.set_rmax(20.0)
ax.grid(True)
ticks = [5,10,15]
ax.set_rticks(ticks)
ax.set_yticklabels(['{0:d} kpc'.format(x) for x in ticks])
ax.set_xlabel("Galactic Longitude", labelpad=15)
ax.tick_params(axis='y', labelsize=20)
# fig.savefig("/Users/adrian/papers/proposals/MDM-2015/GASS.pdf")
# ------
fig,ax = pl.subplots(1,1,figsize=(7,6))
gass_cyl = gass.galactocentric.represent_as(coord.CylindricalRepresentation)
css_cyl = css_c.galactocentric.represent_as(coord.CylindricalRepresentation)
ax.plot(css_cyl.rho.to(u.kpc),
css_cyl.z.to(u.kpc),
color='k', linestyle='none', marker='.', alpha=0.25)
cc = ax.scatter(gass_cyl.rho.to(u.kpc),
gass_cyl.z.to(u.kpc),
c=gass_tbl['VmagAvg'], marker='o')
ax.set_xlim(0,25)
ax.set_ylim(-12.5,12.5)
ax.set_xlabel("R [kpc]")
ax.set_ylabel("z [kpc]")
fig.colorbar(cc)
# ticks = [10,20]
# ax.set_rticks(ticks)
# ax.set_yticklabels(['{0:d} kpc'.format(x) for x in ticks])
# ax.set_xlabel("Galactic Longitude", labelpad=15)
# ax.tick_params(axis='y', labelsize=20)
# # fig.savefig("/Users/adrian/papers/proposals/MDM-2015/GASS.pdf")
In [ ]:
pl.hist(gass_tbl['VmagAvg'])
pl.xlabel("V [mag]")
pl.ylabel("N")
In [ ]:
window_corners = [coord.SkyCoord(l=130.*u.deg, b=20.*u.deg,frame='galactic'),
coord.SkyCoord(l=180*u.deg, b=40*u.deg,frame='galactic')]
css_sky_window_ix2 = coords_in_rect(css_c.galactic, window_corners)
print("{} CSS targets in this window.".format(css_sky_window_ix2.sum()))
In [ ]:
css_distance_ix2 = ((css_c.distance > 11*u.kpc) & (css_c.distance < 33.*u.kpc))
In [ ]:
a13_tbl = css[css_sky_window_ix2 & css_distance_ix2]
a13 = coord.SkyCoord(l=a13_tbl['l']*u.deg, b=a13_tbl['b']*u.deg,
distance=a13_tbl['helio_dist']*u.kpc, frame='galactic')
In [ ]:
print("{} A13 targets".format(len(a13)))
In [ ]:
fig,axes = pl.subplots(1,2,figsize=(12,5))
n,bins,pa = axes[0].hist(a13_tbl['helio_dist'], bins=np.linspace(10,40,12))
axes[0].set_xlabel("Radial dist. [kpc]")
axes[0].set_xlim(bins.min(),bins.max())
# axes[1].hist(gc_cyl_triand.rho, bins=8)
# axes[1].set_xlabel("Cylindrical dist. [kpc]")
# axes[1].set_xlim(17,24)
axes[0].set_ylabel("Number RRL")
In [ ]:
fig,ax = pl.subplots(1,1,figsize=(6,6),subplot_kw =dict(polar=True))
ax.add_artist(mpl.patches.Circle((-8.,0), radius=2.5, transform=ax.transData._b, facecolor='r', alpha=0.2))
ax.add_artist(mpl.patches.Circle((-8.,0), radius=0.5, transform=ax.transData._b, facecolor='y', alpha=1.))
a13_cyl = a13.galactocentric.represent_as(coord.CylindricalRepresentation)
css_cyl = css_c.galactocentric.represent_as(coord.CylindricalRepresentation)
ax.plot(css_cyl.phi.to(u.radian),
css_cyl.rho.to(u.kpc),
color='k', linestyle='none', marker='.', alpha=0.1)
# ax.plot(gass_cyl.phi.to(u.radian), gass_cyl.rho.to(u.kpc),
# color='k', linestyle='none', marker='o')
ax.scatter(a13_cyl.phi.to(u.radian), a13_cyl.rho.to(u.kpc),
c=a13_tbl['VmagAvg'], marker='o')
ax.set_rmax(40.0)
ax.grid(True)
ticks = [10,20,30]
ax.set_rticks(ticks)
ax.set_yticklabels(['{0:d} kpc'.format(x) for x in ticks])
ax.set_xlabel("Galactic Longitude", labelpad=15)
ax.tick_params(axis='y', labelsize=20)
# fig.savefig("/Users/adrian/papers/proposals/MDM-2015/GASS.pdf")
# ------
fig,ax = pl.subplots(1,1,figsize=(7,6))
ax.plot(css_cyl.rho.to(u.kpc),
css_cyl.z.to(u.kpc),
color='k', linestyle='none', marker='.', alpha=0.25)
cc = ax.scatter(a13_cyl.rho.to(u.kpc),
a13_cyl.z.to(u.kpc),
c=a13_tbl['VmagAvg'], marker='o')
ax.set_xlim(0,40)
ax.set_ylim(-20.,20)
ax.set_xlabel("R [kpc]")
ax.set_ylabel("z [kpc]")
fig.colorbar(cc)
# ticks = [10,20]
# ax.set_rticks(ticks)
# ax.set_yticklabels(['{0:d} kpc'.format(x) for x in ticks])
# ax.set_xlabel("Galactic Longitude", labelpad=15)
# ax.tick_params(axis='y', labelsize=20)
# # fig.savefig("/Users/adrian/papers/proposals/MDM-2015/GASS.pdf")
In [ ]:
pl.hist(a13_tbl['VmagAvg'])
pl.xlabel("V [mag]")
pl.ylabel("N")
Overlap between the two samples?
In [ ]:
match_idx,sep,_ = gass.match_to_catalog_sky(a13)
_,sep_a13,_ = a13.match_to_catalog_sky(gass)
In [ ]:
match_idx.shape, gass.shape, sep_a13.shape
In [ ]:
print("{} overlapping targets".format((sep.arcsecond < 1.).sum()))
In [ ]:
from astropy.table import vstack
In [ ]:
gass_tbl['structure'] = 'GASS'
out_gass_tbl = gass_tbl[sep.arcsecond > 1.]
a13_tbl['structure'] = 'A13'
a13_tbl[sep_a13.arcsecond > 1.]['structure'] = 'both'
out_a13_tbl = a13_tbl
In [ ]:
out_tbl = vstack((out_gass_tbl, out_a13_tbl))
out_tbl.sort('ra')
out_tbl['ID'] = ["GA-RR{0:d}".format(x+1) for x in np.arange(len(out_tbl)).astype(int)]
In [ ]:
len(out_tbl)
In [ ]:
ascii.write(out_tbl, "/Users/adrian/projects/triand-rrlyrae/data/targets/GASSA13_2016.txt")
In [ ]:
ascii.write(out_tbl[['ID','ra','dec']], "/Users/adrian/projects/triand-rrlyrae/data/targets/GASSA13_2016_short.txt")#, format="ascii")
In [ ]:
out_tbl_sex = out_tbl.copy()
ra = coord.Longitude(out_tbl_sex['ra']*u.deg)
out_tbl_sex['ra_sex'] = ra.to_string(unit=u.hour, precision=5, sep=':')
dec = coord.Latitude(out_tbl_sex['dec']*u.deg)
out_tbl_sex['dec_sex'] = dec.to_string(unit=u.degree, precision=5, sep=':')
ascii.write(out_tbl_sex[['ID','ra_sex','dec_sex']], "/Users/adrian/projects/triand-rrlyrae/data/targets/GASSA13_2016_short_sexa.txt")
In [ ]:
from gary.observation import distance, apparent_magnitude
from gary.observation.rrlyrae import M_V
In [ ]:
distance(12.5 - M_V(-1.5)).to(u.kpc)
In [ ]:
apparent_magnitude(M_V(-1.5), 10.*u.kpc)
In [ ]:
from scipy.interpolate import InterpolatedUnivariateSpline
In [ ]:
pl.semilogy([14.5,16,17.5], [150,600,2400])
pl.xlim(14, 18)
In [ ]:
def Vmag_to_exptime(V):
s = InterpolatedUnivariateSpline([14.5,16,17.5], np.log10([150,600,2400]), k=1)
y = 10**s(V)
return y
In [ ]:
a13_exptimes = [Vmag_to_exptime(V) for V in a13_tbl['VmagAvg']]
GASS_exptimes = [Vmag_to_exptime(V) for V in gass_tbl['VmagAvg']]
In [ ]:
len(a13_exptimes), len(GASS_exptimes)
In [ ]:
print("A13", (sum(a13_exptimes)*u.second).to(u.hour))
print("GASS", (sum(GASS_exptimes)*u.second).to(u.hour))
In [ ]:
time_grid = astroplan.time_grid_from_range(time_range)
In [ ]:
min_secz = []
min_secz_times = []
for c in a13:
aa = c.transform_to(coord.AltAz(obstime=time_grid, location=mdm.location))
aa = aa[aa.alt > 0.*u.deg]
ix = aa.alt.argmax()
min_secz.append(aa.secz[ix])
min_secz_times.append(time_grid[ix])
for c in gass:
aa = c.transform_to(coord.AltAz(obstime=time_grid, location=mdm.location))
aa = aa[aa.alt > 0.*u.deg]
ix = aa.alt.argmax()
min_secz.append(aa.secz[ix])
min_secz_times.append(time_grid[ix])
In [ ]:
min_secz = u.Quantity(min_secz)
min_secz_times = min_secz_times
In [ ]:
hrs = [mdm.astropy_time_to_datetime(mst).hour for mst in min_secz_times]
In [ ]:
(min_secz < 1.5).sum()
In [ ]:
pl.hist(min_secz)
pl.xlabel("Minimum Airmass")
In [ ]:
pl.hist(hrs)
# pl.xlabel("Minimum Airmass")
In [ ]:
In [ ]: